\section{随机服务系统模拟 *}

我们在 \(§{3.1}\) 中讲到,较为复杂的随机模型往往难以进行彻底的理论分析,这时常常使用随机模拟方法产生模型的大量数据, 从产生的数据对模型进行统计推断。随机服务系统就是这样的一种模型, 经常需要利用随机模拟方法进行研究。

随机服务系统在我们日常生活、工业生产、科学技术、军事领域中是经常遇到的随机模型,比如,研究银行、理发店、商店、海关通道、高速路收费口等服务人员个数的设置和排队规则,研究计算机网络网关、移动网络的调度规则,等等。

{\color{red}\textbf{[实际应用]:} 排队系统无处不在——从银行柜台到网络数据包路由,从医院挂号到CPU任务调度。核心问题都是:如何在有限服务资源下,让顾客等待时间最短、资源利用率最高?}

在概率统计理论中排队论用来研究随机服务系统的数学模型, 可以用来设计适当的服务机构数目和排队规则。如下面的 \(\mathrm{M}/\mathrm{M}/1\) 排队系统。

例 3.4.1. 设某银行仅有一个柜员,并简单假设银行不休息。顾客到来间隔的时间服从独立的指数分布 \(\operatorname{Exp}\left( \lambda \right) (1/\lambda\) 为间隔时间的期望值),如果柜员正在为先前的顾客服务,新到顾客就排队等待,柜员为顾客服务的时间服从均值为 \(1/\mu\) 的指数分布,设 \(u = \lambda /\mu  < 1\) 。设 \({X}_{t}\) 表示 \(t\) 时刻在银行内的顾客数 (包括正在服务的和正在排队的),则 \({X}_{t}\) 是一个连续时马氏链。

{\color{red}\textbf{[解释]:} M/M/1记号表示:第一个M表示到来间隔服从指数分布(M代表无记忆性Markovian),第二个M表示服务时间服从指数分布,1表示只有一个服务台。条件\(u = \lambda/\mu < 1\)确保系统稳定——顾客离开速度快于到来速度。}

这是一个生灭过程马氏链, 有理论结果当系统处于稳定状态时

\[
P\left( {{X}_{t} = i}\right)  = {u}^{i}\left( {1 - u}\right) ,i = 0,1,2,\ldots
\]

设随机变量 \(N\) 服从 \({X}_{t}\) 的平稳分布。于是银行中平均顾客数为

\[
{EN} = \frac{u}{1 - u}
\]

{\color{red}\textbf{[物理图像]:} 当\(u = 0.9\)时,\(EN = 9\)——虽然柜员只忙90\%的时间,但平均有9人在银行!这反映了排队系统的``拥堵放大效应'':利用率越接近100\%,等待队列呈指数增长。}

平均队列长度 \({EQ}\) 等于 \({EN}\) 减去平均正在服务人数,正在服务人数 \({Y}_{t}\) 为

\[
{Y}_{t} = \left\{  \begin{array}{ll} 1, & \text{ 当 }{X}_{t} > 0, \\  0, & \text{ 当 }{X}_{t} = 0 \end{array}\right.
\]

所以 \(E{Y}_{t} = P\left( {{Y}_{t} = 1}\right)  = 1 - P\left( {N = 0}\right)  = u\) ,于是平均队列长度为

\[
{EQ} = {EN} - E{Y}_{t} = \frac{u}{1 - u} - u = \frac{{u}^{2}}{1 - u},
\]

设顾客平均滞留时间为 \({ER}\) ,由关系式

\[
{EN} = \lambda  \cdot  {ER}
\]

{\color{red}\textbf{[关键理解]:} 这是Little定律——排队论的基本公式:平均人数=到达率×平均停留时间。直观理解:如果每秒到来\(\lambda\)个顾客,每人平均停留\(ER\)秒,那么任意时刻平均有\(\lambda \cdot ER\)人在系统中。}

可知平均滞留时间为

\[
{ER} = \frac{u}{\lambda \left( {1 - u}\right) } = \frac{1}{\mu  - \lambda },
\]

进一步分析还可以知道顾客滞留时间 \(R\) 服从均值为 \(1/\left( {\mu  - \lambda }\right)\) 的指数分布。

从上面的例子可以发现, 一个随机服务系统的模型应该包括如下三个要素:

\begin{itemize}
\item 输入过程: 比如, 银行的顾客到来的规律。
\end{itemize}

\begin{itemize}
\item 排队规则: 比如, 银行有多个柜员时, 顾客是选最短一队, 还是随机选一队, 还是统一排队等候叫号, 顾客等得时间过长后是否会以一定概率放弃排队。
\end{itemize}

\begin{itemize}
\item 服务机构:有多少个柜员,服务时间的分布等。
\end{itemize}

虽然某些随机服务系统可以进行严格理论分析得到各种问题的理论解, 但是, 随机服务系统中存在大量随机因素, 使得理论分析变得很困难以至于不可能。例如, 即使是上面的银行服务问题, 可能的变化因素就包括: 顾客到来用齐次泊松过程还是非齐次泊松过程, 柜员有多少个, 是否不同时间段柜员个数有变化, 柜员服务时间服从什么样的分布, 顾客排队按照什么规则,是否 VIP 顾客提前服务,顾客等候过长时会不会放弃排队,等等。包含了这么多复杂因素的随机服务系统的理论分析会变得异常复杂, 完全靠理论分析无法解决问题, 这时, 可以用随机模拟方法给出答案。

{\color{red}\textbf{[为什么需要模拟]:} 现实中的排队系统极其复杂——早高峰时银行顾客到来加速,VIP客户插队,顾客等太久会离开。这些因素让解析公式失效,但模拟可以轻松处理:只需在代码中加入相应规则即可。}

在模拟随机服务系统时, 可以按时间顺序记录发生的事件, 如顾客到来、顾客接受服务、 顾客结束服务等,这样的系统的模拟也叫做离散事件模拟。

离散事件模拟算法可以分为三类:

\begin{itemize}
\item 活动模拟, 把连续时间离散化为很小的单位, 比如平均几秒发生一个新事件时可以把时间精确到毫秒, 然后时钟每隔一个小时间单位前进一步并查看所有活动并据此更新系统状态。缺点是速度太慢。
\end{itemize}

\begin{itemize}
\item 事件模拟。仅在事件发生时更新时钟和待发生的事件集合。这样的方法不受编程语言功能的限制, 运行速度很快, 也比较灵活, 可以实现复杂逻辑, 但是需要自己管理的数据结构和逻辑结构比较复杂,算法编制相对较难。
\end{itemize}

\begin{itemize}
\item 过程模拟。把将要到来的各种事件看成是不同的过程, 让不同的过程并行地发生, 过程之间可以交换消息, 并在特殊的软件包或编程语言支持下自动更新系统状态。在计算机实现中需要借助于线程或与线程相似的程序功能。这类软件包有 \(\mathrm{C} +  +\) 语言软件包 \(\mathrm{C} +  +\) SIM 和 Python 语言软件包 SimPy。优点是系统逻辑的编码很直观,程序模块化, 需要用户自己管理的数据结构和逻辑结构少。过程模拟是现在更受欢迎的离散事件模拟方式。
\end{itemize}

{\color{red}\textbf{[类比]:} 三种方法类比观察交通:活动模拟像每秒拍一张照片(大量无用帧);事件模拟像只在车辆通过路口时记录(高效但需手动追踪所有车);过程模拟像给每辆车装GPS自动上传(最直观,软件自动管理)。}

事件模拟算法必须考虑的变量包括:

\begin{itemize}
\item 当前时刻 \(t\) ;
\end{itemize}

\begin{itemize}
\item 随时间而变化的计数变量,如 \(t\) 时刻时到来顾客人数、已离开人数;
\end{itemize}

\begin{itemize}
\item 系统状态, 比如是否有顾客正在接受服务、排队人数、队列中顾客序号。
\end{itemize}

这些变量仅在有事件发生时 (如顾客到来、顾客离开) 才需要记录并更新。为了持续模拟后续事件, 需要维护一个将要发生的事件的列表 (下一个到来时刻、下一个离开时刻), 列表仅需要在事件发生时进行更新。其它变量可以从这三种变量中推算出来,比如,顾客 \(i\) 在时间 \({t}_{1}\) 时到达并排队,在时间 \({t}_{2}\) 时开始接受服务,在时间 \({t}_{4}\) 时结束服务离开,则顾客 \(i\) 的滞留时间为 \({t}_{4} - {t}_{1}\) 。

例 3.4.2. 用事件模拟的方法来模拟例3.4.1。目的是估计平均滞留时间 \({ER}\) 。想法是,模拟生成各种事件发生的时间,模拟很长时间,丢弃开始的一段时间 \({T}_{0}\) 后,用 \({T}_{0}\) 后到达的顾客的总滞留时间除以 \({T}_{0}\) 后到达的顾客人数来估计平均滞留时间。设 \({T}_{0}\) 时间后每位顾客滞留时间的模拟值为 \({R}_{i},i = 1,2,\ldots ,m\) ,可以用 \(\left\{  {R}_{i}\right\}\) 作为随机变量 \(R\) 的样本来检验 \(R\) 的分布是否指数分布。

{\color{red}\textbf{[关键技巧]:} 丢弃初始时段\(T_0\)称为``预热''(warm-up)——系统从空状态启动需要时间才能进入稳态。就像开车上高速,前几分钟在加速,不能用来估计巡航速度。只有稳态后的数据才符合理论分布。}

用事件模拟方法进行离散事件模拟的算法关键在于计算系统状态改变的时间,即各个事件的发生时间, 这个例子中就是顾客到来、顾客开始接受服务、顾客离开这样三种事件, 由此还可以得到每个顾客排队的时间和服务的时间。

在没有明确算法构思时, 可以从时间 0 开始在纸上按照模型规定人为地生成一些事件并人为地找到下一事件。这样可以找到要更新的数据结构和更新的程序逻辑。

下面的算法保持了一个将要发生的事件的集合, 在每次事件发生时更新时钟, 更新时钟时从事件集合中找到最早发生的事件进行处理, 并生成下一事件到事件集合中, 如此重复直到需要模拟的时间长度。

\HRule

\{初始化当前时钟 \(t \leftarrow  0\) ,柜员忙标志 \(B \leftarrow  0\) ,当前排队人数 \(L \leftarrow  0\) ,

最新到来顾客序号 \(i \leftarrow  0\) ,正在服务顾客序号 \(j \leftarrow  0\) ,已服务顾客数 \(n \leftarrow  0\}\)

从 \(\operatorname{Exp}\left( \lambda \right)\) 抽取 \(X\) ,设置下一顾客来到时间 \(A \leftarrow  X\)

repeat \{

\hspace*{1em} if \(\left( {B = 0\text{ or }\left( {B = 1\text{ and }A < E}\right) }\right) \left\{  {\# E}\right.\) 是正在服务的顾客结束时刻

\hspace*{2em} \(t \leftarrow  A\)

\hspace*{1em} \} else \{

\hspace*{2em} \(t \leftarrow  E\)

\hspace*{1em} \}

\hspace*{1em} if \(\left( {t > {T}_{1}}\right)\) break \(\# {T}_{1}\) 是预先确定的模拟时长

\hspace*{1em} if \(\left( {t =  = A}\right) \{ \#\) 待处理到达事件

\hspace*{2em} \(L \leftarrow  L + 1\)

\hspace*{2em} \(i \leftarrow  i + 1\) ,记录第 \(i\) 位顾客到来时间 \({a}_{i} \leftarrow  t\)

\hspace*{2em} 从 \(\operatorname{Exp}\left( \lambda \right)\) 抽取 \(X,A \leftarrow  t + X\)

\hspace*{2em} if \(\left( {B = 0}\right) \{ \#\) 不用排队,直接服务

\hspace*{3em} \(B \leftarrow  1,L \leftarrow  L - 1\)

\hspace*{3em} \(j \leftarrow  j + 1\) ,置第 \(j\) 位顾客开始服务时间 \({s}_{j} \leftarrow  t\)

\hspace*{3em} 从 \(\operatorname{Exp}\left( \mu \right)\) 抽取 \(Y\) ,置 \(E \leftarrow  t + Y\)

\hspace*{2em} \}

\hspace*{1em} \} else \{ \# 待处理结束服务事件

\hspace*{2em} \(B \leftarrow  0\)

\hspace*{2em} \(n \leftarrow  n + 1\) ,记录第 \(n\) 个顾客结束服务时间 \({e}_{n} \leftarrow  t\)

\hspace*{2em} if \(\left( {L > 0}\right) \{ \#\) 排队顾客开始服务

\hspace*{3em} \(L \leftarrow  L - 1\)

\hspace*{3em} \(B \leftarrow  1\)

\hspace*{3em} \(j \leftarrow  j + 1,{s}_{j} \leftarrow  t\)

\hspace*{2em} 从 \(\operatorname{Exp}\left( \mu \right)\) 抽取 \(Y\) ,置 \(E \leftarrow  t + Y\)

\hspace*{1em} \}

\}

\HRule

\}

\{令 \(I = \left\{  {i : {T}_{0} \leq  {s}_{i} \leq  {T}_{1}}\right\}\) ,求 \(\left\{  {{e}_{i} - {a}_{i},i \in  I}\right\}\) 的平均值作为 \({ER}\) 估计 \(\}\)

{\color{red}\textbf{[解释]:} 事件模拟的核心思想是``时钟跳跃''——时间不连续流动,而是直接跳到下一个有意义的时刻。就像看录像只看关键镜头,跳过无事发生的部分,大幅提高效率。}

从这个算法可以看出, 事件模拟方法需要用户自己管理待处理事件集合与时钟, 算法设计难度较大。

用随机服务系统进行建模和模拟研究的一般步骤如下:

\begin{itemize}
\item 提出问题。比如, 银行中顾客平均滞留时间与相关参数的关系。
\end{itemize}

\begin{itemize}
\item 建立模型。比如, 设顾客到来服从泊松过程, 设每个顾客服务时间服从独立的指数分布。
\end{itemize}

\begin{itemize}
\item 数据收集与处理。比如,估计银行顾客到来速率,每个顾客平均服务时间,等等。
\end{itemize}

\begin{itemize}
\item 建立模拟程序。一般使用专用的模拟软件包或专用模拟语言编程。
\end{itemize}

\begin{itemize}
\item 模拟模型的正确性确认。
\end{itemize}

\begin{itemize}
\item 模拟试验和模拟结果分析。
\end{itemize}

离散事件模拟问题一般都比较复杂, 即使借助于专用模拟软件或软件包, 也很难确保实现的模拟算法与问题的实际设定是一致的。为此, 需要遵循一些提高算法可信度的一般规则。 算法程序一定要仔细检查, 避免出现参数错误、逻辑错误。在开始阶段, 可以利用较详尽的输出跟踪模拟运行一段时间,人工检查系统运行符合原始设定。尽可能利用模块化程序设计, 比如,在 \(\mathrm{M}/\mathrm{M}/1\) 问题模拟中,顾客到来可能遵循不同的规律,比如时齐泊松过程、非时齐泊松过程, 把产生顾客到来时刻的程序片段模块化并单独检查验证, 就可以避免在这部分出错。问题实际设定可能比较复杂, 在程序模块化以后, 如果有一种较简单的设定可以得到理论结果, 就可以用理论结果验证算法输出, 保证程序框架正确后, 再利用模块化设计修改各个模块适应实际复杂设定。
